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TENSORIAL DIFFUSION EQUATION 



K. B. NAKSHATRALA AND A. J. VALOCCHI 

Abstract. We consider the tensorial diffusion equation, and address tiie discrete maximum- 
minimum principie of mixed finite eiement formulations. In particular, we address non-negative 
solutions (wfiidi is a special case of tfie maximum-minimum principle) of mixed finite element formu- 
lations. It is well-known that the classical finite element formulations (like the single-field Galerkin 
formulation, and Raviart- Thomas, variational multiscale, and Galerkin/least-squares mixed formu- 
lations) do not produce non-negative solutions (that is, they do not satisfy the discrete maximum- 
minimum principle) on arbitrary meshes and for strongly anisotropic diffusivity coefficients. 

In this paper we present two non-negative mixed finite element formulations for tensorial diffu- 
sion equations based on constrained optimization techniques. These proposed mixed formulations 
produce non-negative numerical solutions on arbitrary meshes for low-order (i.e., linear, bilinear 
and trilinear) finite elements. The first formulation is based on the Raviart-Thomas spaces, and 
the second non-negative formulation is based on the variational multiscale formulation. For the 
former formulation we comment on the effect of adding the non-negative constraint on the local 
mass balance property of the Raviart-Thomas formulation. 

We perform numerical convergence analysis of the proposed optimization-based non-negative 
mixed formulations. We also study the performance of the active set strategy for solving the 
resulting constrained optimization problems. The overall performance of the proposed formulation 
is illustrated on three canonical test problems. 

1. INTRODUCTION 

Robustness of numerical methods for flow and transport problems in porous media is important 
for development of simulators to be used in a wide range of applications in subsurface hydrology and 
contaminant transport. In order to obtain robust and reliable numerical results it is imperative 
to preserve basic properties of solutions of mathematical models by computed approximations. 
In the simulation of reactive transport of contaminants one such basic properties is non-negative 
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solutions as concentration of a chemical or biological species physically can never be negative. Since 
the domain of interest in subsurface flows is highly complex, one needs to employ unstructured 
computational grids. Therefore, obtaining non-negative solutions on unstructured meshes is an 
essential feature in the simulation of reactive transport of contaminants, as well as many other 
physical processes. 

However, obtaining non- negative solutions on unstructured grids using a numerical method (fi- 
nite element, finite volume or finite difference) is not an easy task. In addition, there is another 
complexity arising from an anisotropic diffusion tensor. In subsurface flows, the heterogeneity in 
the velocity field will give rise to a non-homogeneous anisotropic diffusion tensor with non-negligible 
cross terms [1]. Several studies have shown that standard treatment of the cross-diffusion term will 
result in negative solutions on general computational grids, see |2j and references therein. Several 
ad-hoc procedures have been proposed in the literature. For example, a post processing step is 
typically employed in which one performs some sort of "smoothing." But this procedure, in many 
cases, is not variationally consistent. Some other methods are limited in their range of applicability 
(e.g., the method proposed in Reference [2] can handle only structured grids). 

Herein we consider the dispersion/diffusion process for steady single-phase flow in heterogeneous 
anisotropic porous media. Such a flow can be described by the Poisson's equation with a tensorial 
diffusion coefficient, which when written in the mixed form is similar to the governing equations of 
Darcy flow. In this paper we propose two optimization-based mixed methods for solving tensorial 
diffusion equation that gives non-negative solutions on general grids for linear finite elements. The 
two methods are developed by rewriting the Raviart-Thomas and variational multiscale formula- 
tions as constrained minimization problems subject to a constraint on the primary variable to be 
non-negative. A similar approach based on optimization techniques has been used by Liska and 
Shashkov [3J for a single field formulation, but herein we consider mixed finite element formulations. 

The main idea behind the proposed methods is to augment a constraint on nodal values to be non- 
negative to the discrete (that is, after spatial finite element discretization) variational statement of 
the underlying formulation. Since we consider only low- order finite elements (and since the shape 
functions for these elements do not change their sign within an element), non-negative nodal values 
ensure non-negative solution every where in the element, and hence non-negative solution on the 
whole domain. This argument will not hold for high-order finite elements as shape functions for 
these finite elements (in general) change sign within an element. 
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Throughout this paper continuum vectors are denoted with lower case boldface normal letters, 
and (continuum) second-order tensors will be denoted using (ETeX) blackboard font (for example, 
V and D, respectively). We denote finite element vectors and matrices with lower and upper case 
boldface italic letters, respectively. For example, vector v and matrix K. The curled inequality 
symbols ^ and ^ are used to denote generalized inequalities between vectors, which represent 
component- wise inequalities. That is, given two vectors a and b, a y b means Oj > bi A 
similar definition holds for the symbol ^. Other notational conventions adopted in this paper are 
introduced as needed. 

1.1. Governing equations. Let C R'"^ (where "nd" is the number of spatial dimensions) be 
a bounded domain with boundary dVL = \ 0, where denotes the closure of J7. Consider the 
diffusion of a chemical species in anisotropic heterogeneous medium, which is governed by the 
second-order elliptic tensorial diffusion partial differential equation. The governing equations are 

(1) -V • (]D)(x)Vc(x)) = /(x) in 

(2) -n(x) • D(x)Vc = tP(x) on 

(3) c(x)=cP(x) onr° 

where c(x) denotes the concentration field, /(x) is the volumetric source, tP(x) is the prescribed flux 
(i.e., Neumann boundary condition), cP(x) is the prescribed concentration (i.e., Dirichlet boundary 
condition), T^ is that part of the boundary on which Dirichlet boundary condition is applied, 
is the part of the boundary on which Neumann boundary condition is applied, n(x) is unit 
outward normal to the boundary, and V denotes gradient operator. For well-posedness one requires 
-pD y -pN _ .^^^ pD ^ pN _ 0^ g^j^^ £qj, uniqueness F^ ^ 0. We assume that the coefficient of 
diffusivity D(x) is a symmetric positive definite tensor such that, for some < ai < Q2 < +oo, we 
have 

(4) aiy^y < y^B(x)y < aay^y Vx G f], Vy / G M'^'^ 
In addition, we assume that B(x) is continuously differentiable. 

1.2. First-order (or mixed) form. In many situations, the primary quantity of interest is the 
flux. But a single field (or primal) formulation does not produce accurate solutions for the flux. 
One can calculate the flux by differentiating the obtained c(x), but there will be a loss of accuracy 
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during this process. For example, under a single field formulation, linear finite elements produce 
fluxes that are constant and discontinuous across elements. This means that there is no flux balance 
across element edges. Balance of flux along element edges is a highly desirable feature and is of 
physical importance in many practical engineering problems. In order to alleviate aforementioned 
drawbacks of single formulations, mixed formulations are often employed. Equations ([I|)-(l3|) in 
mixed (or first-order) form can be written as 

(5) D~^(x)v(x) = -Vc inn 

(6) V • V = /(x) in Q 

(7) v(x) • n(x) = tP(x) onT^ 

(8) c(x) = cP(x) on r° 

where v(x) is an auxiliary variable, which can be interpreted as follows: given a plane defined by 
a normal n, the quantity v • n will be the flux through the plane. 



1.3. Maximum-minimum principle. It is well-known that some elliptic partial differential equa- 
tions (under appropriate regularity assumptions) satisfy the so-called maximum-minimum principle, 
and the Poisson's equation is one of them [3]. We now state the classical maximum- minimum prin- 
ciple for second-order elliptic partial differential equations. (In Section [3] we state and prove a 
maximum- minimum principle under milder regularity assumptions. Note that weak solutions may 
also possess a maximum- minimum principle. For example, see Reference [5].) Consider the bound- 
ary value problem given by equations ([I])-®. Let c(x) G C^(r2) n C^{Q), where C^(0) denotes the 
set of twice continuously differentiable functions defined on J7, and C^{Cl) the set of uniformly con- 
tinuous functions defined on Q. If /(x) > (or if /(x) < 0) in 17 then c(x) attains its minimum (or 
its maximum) on the boundary of 17. For a detailed discussion on maximum-minimum principles 
see References [H El O [5] . 

One of the important consequences of maximum-minimum principles is the non-negative solution 
of a (tensorial) diffusion equation under non-negative forcing function with non-negative prescribed 
Dirichlet boundary condition. Obtaining non-negative solutions is of paramount importance in 
studying transport of chemical and biological species as negative concentration of a species is 
unphysical. 

4 



1.4. Discrete maximum-minimum principle. The discrete analogy of the maximum-minimum 
principle is commonly referred to as the discrete maximum-minimum principle (DMP). However, 
many numerical formulations do not unconditionally satisfy the discrete maximum-minimum prin- 
ciple. Typically, there will be restrictions on the mesh or on the magnitude of coefficients of the 
diffusivity tensor. For example, the single field Galerkin formulation in the case of scalar diffusion 
satisfies the discrete maximum-minimum principle if the mesh satisfies weak acute condition [8] 
(and also see Appendix). The question whether we get non- negative numerical solutions leads us 
to the discrete maximum-minimum principle. 

For recent works on DMP see References [H [lOl [TTl [121 El EH US] and also see the discussion in 
Reference [3l Introduction]. Considerable attention to DMP has also been given in the finite volume 
literature [21 [HI [171 [18] . Optimization-based techniques have been employed in References and 
[19| to address DMP. For completeness, some of the classical results on discrete maximum-minimum 
principle are outlined in Appendix. 

In this paper we concentrate on obtaining non-negative numerical solutions using mixed formu- 
lations under non-negative forcing function with non-negative Dirichlet boundary condition where 
ever it is prescribed. (Note that we do not assume that the Dirichlet boundary condition has to be 
prescribed on the whole boundary.) In all our test problems (see Section [4]) we have /(x) > in $7 
and cP(x) > on dO,. By using the maximum-minimum principle one can conclude that c(x) > 
in whole of ft (the closure of Q). That is, for the chosen test problems in Section [H we must have 
non-negative solutions in the whole domain. 

1.5. Main contributions of this paper. Some of the main contributions of this paper are as 
follows: 

• We numerically demonstrate that various conditions outlined in Appendix (which are suf- 
ficient for isotropic diffusion) are not sufficient for tensorial diffusion equation under the 
Raviart-Thomas and variational multiscale formulations to produce non-negative solutions 
under non-negative forcing functions with non-negative prescribed Dirichlet boundary con- 
ditions. 

• We develop a non-negative formulation based on the lowest-order Raviart-Thomas spaces, 
and discuss the consequences of obtaining non-negative solutions on the local mass balance 
property. 
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• We extend the variational multiscale formulation to produce non-negative solutions on gen- 
eral grids for low-order finite elements under non-negative forcing function and prescribed 
non-negative Dirichlet boundary condition. We also show that the (continuous) variational 
multiscale formulation satisfies a continuous maximum-minimum principle (under appro- 
priate regularity assumptions). 

1.6. Organization of the paper. The remainder of this paper is organized as follows. In Section 
[2] we present a non-negative formulation based on the lowest-order Raviart-Thomas (RTO) finite 
element spaces, which is achieved by adding a non- negative constraint to the discrete variational 
setting of the Raviart-Thomas formulation. For this non-negative formulation we present both 
primal and dual constrained optimization problems, and comment on the ease of solving these 
problems and also the consequences of imposing the non-negative constraint on the local mass 
balance. In Section [3l a non-negative formulation based on the variational multiscale formulation 
will be presented. We also show that the (continuous) variational multiscale formulation satisfies 
a continuous maximum-minimum principle. Numerical results along with a discussion on the nu- 
merical performance of both the proposed non-negative formulations will be presented in Section 
m Finally, conclusions are drawn in Section O 



2. A NON-NEGATIVE MIXED FORMULATION BASED ON RAVIART-THOMAS 

SPACES 

The Raviart-Thomas finite element formulation is widely used (for an example in subsurface 
modeling see [20]) to solve diffusion equations in mixed form, and is based on the classical mixed 
formulation [21]. The simplest and lowest order Raviart-Thomas space (commonly denoted as RTO) 
consists of fluxes evaluated on the midpoints of edges and constant pressure over elements. We first 
present the weak form and variational structure behind the Raviart-Thomas formulation. We then 
modify the variational structure by adding non-negative constraint on the concentration to build a 
non- negative low-order finite element formulation based on the RTO spaces. 
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To this end, define function spaces as 



(9) 


V 


= {v(x) 


v(x) G 


V-vGL2(J1), v(x) 


• n(x) = tP(x) on 


(10) 


w 


= {v(x) 


v(x) G 


V- V G L^{a), v(x) 


• n(x) = on 


(11) 


V 


= L^{n) 









Recall that "nd" denotes the number of spatial dimensions. For further details on function spaces 
see the monograph by Brezzi and Fortin [22]. Let w(x) and (/(x) denote the weighting functions 
corresponding to v(x) and c(x), respectively. The classical mixed formulation for equations (IS])-® 
can be written as: Find v(x) G V and c(x) G V such that 

(12) (w;D"V) -(V-w;c) + (w-n;cP)pD-(g;V-v-/)=0 Vw(x) G W, g(x) G P 

It is well-known that, under appropriate smoothness conditions on the domain and its boundary, 
the above saddle-point formulation is well-posed [22] • That is, a unique (weak) solution exists 
for this problem that depends continuously on the input data. However, to obtain stable results 
using a finite element approximation, the finite dimensional spaces dV and d V in which 
a numerical solution is sought have to satisfy the Ladyzhenskaya-Babuska-Brezzi (LBB) stability 
condition [22] . One such space that satisfies the LBB condition is the popular Raviart-Thomas (RT) 
finite element space. In this paper we consider only the lowest order Raviart-Thomas triangular 
finite element space (RTO). Let 7/j be a triangulation on $7. The lowest order Raviart-Thomas 
finite dimensional subspaces on triangles are defined as 

(13) Vh '■= {p \ P = 9- constant on each triangle K G T^} 

(14) Vh := {v = (t;«,i;(2)) | v'^^^ = ax + bKX, = ck + bKy;aK,bK, ck e R; K £ 

2.1. Discrete equations. The discretized finite element equations of the Raviart-Thomas formu- 
lation for the mixed form of tensorial diffusion equation can be written as [20] 



(15) 


VV 




1 




hi 


hA 




^pv 


O 











where O is a zero matrix of appropriate size, the matrix is symmetric and positive defi- 
nite, V denotes the (finite element) vector of flux degrees-of-freedom, and p denotes the vector of 
concentration degrees-of-freedom. Comparing the weak form (I12p and discrete equation (llSp . the 
matrices Km, and Kpy are obtained after the finite element discretization of the terms (w; D~"'^v) 
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and — (g; V • v), respectively. Since equation (jl2p is written in symmetric form, K^tp (which comes 
from the term — (V • w; c)) will be equal to K^^. The vectors and fp are, respectively, obtained 
from — (w • n;cP)pD and —{q;f) after the finite element discretization. Since there is no term in 
equation (jl2p that contains both q in the weighting (i.e., first) slot and c in the second slot in a 
bilinear form {■;■), we have the matrix Kpp = O (a zero matrix). 

The above system of equations (|15|) is equivalent to the following constrained minimization 
problem 



(16) (Pl-RTO) 



minimize^ ^v'^Ki,i,v — v'^ 
subject to KpyV — fp = 

where is a zero vector of appropriate size. Note that the constraint in equation (I16p is the local 
mass balance condition for each element. We refer the above equation as the primal problem for 
the Raviart-Thomas formulation, and denote it as (Pl-RTO). This primal problem belongs to the 
class of convex quadratic programming problems, and from optimization theory (for example, see 
Reference [23]) it can be shown that the problem has a unique global minimizer. 

Remark 2.1. A quadratic program is an optimization problem in which the objective function is a 
quadratic function and the (equality and inequality) constraints are all linear. In a convex quadratic 
program the Hessian of the objective function is positive semidefinite. 

Remark 2.2. It is interesting to note that, from the complexity theory, problem (jl6p can be solved 
in polynomial time (for example, using the ellipsoid and interior point methods) [231 I24j . Note 
that the term "polynomial time" in the context of complexity theory should not be confused with the 
term "polynomial convergence, " which is commonly used in the convergence studies using the finite 
element method. 

Define the Lagrangian as 
(17) C{v,p):= ^v^K,,v - v^f„ + [Kp,v - fp) 

where p is the vector of Lagrange multipliers. Using the Lagrange multiplier method [23] the primal 
problem (jl6p is equivalent to 



(18) extremize C{v,p) 



and the first-order optimality conditions for tliis problem gives rise to the discretized finite element 
equations (jlSp . We now write the dual problem corresponding to the primal problem ()16p . To this 
end, define the Lagrange dual function as 

(19) g{p) := inf C[v,p) = -^p^Kp^R-^^K^, p + p^ {K,,K-^f, - f^) - \flK-^f, 

The above expression on the right-hand side is obtained as follows. Let v* be the minimizer that 
gives the infimum of C{v,p) with respect to v. Then v* has to satisfy 

(20) K^^v* -f, + Kl^p = 

which is a necessary condition, and is obtained by equating the derivative of C{v,p) (which is 
defined in equation (fT7|) ) with respect to v to zero. Since the matrix K^y is positive definite (and 
hence invertible) we have 

(21) V* = (/, - K» 

By substituting the above expression for the minimizer v* into the definition of C{v,p) (|17p we 
obtain the expression on the right-hand side of equation (jl9p . 

The dual problem corresponding to the primal problem (jl6p can then be written as 

(22) maximize g{p) 



p 



which is equivalent to 

(23) (Dl-RTO) minimize \p^ Kp^R-^ K^^ p - p^ {Kp^R-^ f , - f ;) 
The stationarity of the above problem implies 

(24) KpyKyy -^pv P ~ -^pv^vv fv ~ fp 

which is the Schur complement form of equation (llSp expressed in terms of Lagrange multipliers 

pv-'-^vv ^^pv 



by analytically eliminating the variable v. Note that the Schur complement operator Kp^K^^K"^ 



is symmetric and positive definite. 

A simple numerical example to be presented later (e.g., see Figure ED shows that RTO triangular 
element does not satisfy the discrete maximum-minimum principle, and in particular, does not 
produce non-negative solutions for non-negative forcing functions with non-negative prescribed 
Dirichlet boundary conditions. 



2.2. A non-negative mixed formulation. In order to get non-negative solutions under the RTO 
spaces, we pose the dual problem as 



(25) (D2-RT0) 



minimizcp ^p'^ Kp^KjK'^^ P - {Kp^Kjf^ - f^) 
subject to p h 



Recall that the symbol >z denotes the generalized inequality between vectors, which represents 
component-wise inequality (see Introduction, just above subsection ll.H for a discussion on this 
notation). The primal problem corresponding to this new dual problem will then be 



minimize^ ^v'^KyyV — v'^ 
subject to Kpi,v — fp ^ 



(26) (P2-RT0) 



Remark 2.3. The primal problem given in equation (j26p is obtained by inspection. That is, one 
can easily check (using a direct calculation) that the dual problem of this new primal problem ()26p 
will be the same as equation ()25p . Also, it should be noted that one can write the dual problem 
corresponding to a given dual problem (that is, the dual of a dual). For the problem at hand, the 
dual of the dual problem will be the same as the primal problem, which is not the case in general 
|23j . Hence, one will obtain the primal problem ()26p by writing the dual of the dual problem (j25p . 

By comparing the constraints in equations (jl6p and (|26p one can conclude that under the pro- 
posed non-negative method based on the Raviart-Thomas formulation one may violate local mass 
balance by creating artificial sinks. We can infer more on local mass balance by looking at the 
Karush-Kuhn- Tucker (KKT) conditions (which in this case are necessary and sufficient for the op- 
timality) for the new primal problem given by equations (j26p . The KKT optimality conditions for 
the (P2-RT0) problem are 

(27) K,,v + Kl,p = f, 

(28) Kp^v -fp^O 

(29) ptO 

(30) piiKp^v- fp)i = Vi 

The last condition (which is basically the complementary slackness condition in the KKT system 

of equations) implies that one may not have local mass balance in those elements for which the 
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Lagrange multiplier vanishes (i.e., pi = 0, where i denotes the element number). Note that in the 
RTO formulation, the Lagrange multiplier pi denotes the concentration in the i^^ element. 

Remark 2.4. From optimization theory [53] one can show that the primal (P2-RT0) and dual 
(D2-RT0) problems are equivalent. That is, there is no duality gap for the optimization-based RTO 
formulation. The difference between primal and dual solutions is commonly referred to as the duality 
gap. In general, the solution of a dual problem gives an upper bound to its corresponding primal 
problem. 

However, from a computational point of view the primal and dual problems can have different 
numerical performance (especially with respect to computational cost, and selection of numerical 
solvers). The primal problem (P2-RT0) has more complicated constraints than the dual problem 
(D2-RT0) for which the constraints are (lower) bounds on the design variable p. Compared to the 
primal problem (P2-RT0), the dual problem (D2-RT0) has a more complicated objective function, 
which is defined in terms of Schur complement operator. For all the numerical results presented in 
this paper, we have used the dual problem (D2-RT0). However, we have compared the numerical 
solutions obtained using the dual problem with the primal problem (P2-RT0), and the solutions are 
identical as predicted by the theory. 

Special solvers are available in the literature (for example, especially designed interior point meth- 
ods [25^124] ] that are effective for solving problems that belong to the class of quadratic programming 
with constraints being just bounds on design variables. Similarly, special solvers are available for 
problems involving Schur complement operators. For example, the preconditioned conjugate gradient 
(PCG) solver is quite effective for solving large-scale problems involving Schur complement opera- 
tor. A detailed analysis comparing computational costs of the primal (P2-RT0) and dual (Q2-RT0) 
problems is beyond the scope of this paper. 

3. A NON-NEGATIVE VARIATIONAL MULTISCALE MIXED FORMULATION 

Masud and Hughes [26] have proposed a stabilized mixed formulation for the first-order form of 

the Poisson equation that satisfies the LBB condition. In this paper we refer to this formulation 

as the variational multiscale (VMS) formulation. Nakshatrala et al. [2^ have shown that the 

variational multiscale formulation can be derived based on the multiscale framework proposed by 

Hughes [28]. The variational multiscale formulation possesses many favorable numerical properties 

and performs very well in practice. For example, the formulation passes three dimensional patch 
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tests even for distorted elements [27] . Another feature of this formulation worth mentioning is that 
the equal-order interpolation for c and v is stable \28\ I27j . However, the variational multiscale 
formulation in general does not satisfy the discrete maximum-minimum principle, which will be 
illustrated below in Figure [12] using a simple numerical example. In this section, we present a 
non-negative method based on the variational multiscale mixed formulation. To this end, we first 
present the weak form and variational structure behind the variational multiscale formulation. 

Let the function spaces for the concentration and its corresponding weighting function under the 
VMS formulation be 

(31) V = H^{9) 

(32) Q = H^{Q) 

where H^{^) is a standard Sobolev space defined on domain [22]. The variational multiscale 
formulation reads [271 (23: Find c(x) G V and v(x) G V such that 

(w; D-V) - (V • w; c) + (w • n; cP)pD - (<?; V • v - /) 

(33) - ^ (]D)"V + Vg;]D)(]D)"V + Vc)) = Vg(x) G Q, w(x) G W 

where w(x) and (/(x) are weighting functions corresponding to v(x) and c(x), and V and W are 
defined in equations (|9]) and (jlOp . respectively. The stationarity (minimizing with respect to v and 
maximizing with respect to c) of the following (continuous) optimization problem 

(34) extremize - (v; D~ V) - (c; V • v - /) + (v n; cP)rD - - (D~ V + Vc; D(D~ V + Vc)) 

vgv, c&v 2 4 

is equivalent to the weak form given by equation ([33]) . 

Many practically important problems do not have solutions in C^(r2)nC'^(J^), and hence for these 
problem one cannot employ the classical maximum-minimum principle, which we have outlined in 
Section [1] For example, there exist no (classical) solutions to test problems ^1 and 7^2 (which 
are defined in Section [4]) that belong to 0^(0) as the forcing functions in both these cases are not 
continuous on Q.. (On the other hand, the solution to test problem ^^3 does belong to C^(r2) n 
C^{Cl).) However, weak solutions do exist for test problems #1 and #2. Using regularity 
theory (for example, see Reference [29]), one can show that these solutions, in fact, belong to 
H^{Vt)r^C^{^^)f^C^{^l). 
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Remark 3.1. Note that H'^{Q) is not a subset ofC^{^) or vice-versa. On the other hand, H'^{Q) C 

3.1. Continuous maximum-minimum principle. In this subsection we demonstrate one of the 
main contributions of this paper, namely that the weak solution under the variational multiscale 
formulation (under appropriate regularity assumptions) satisfies a continuous maximum-minimum 
principle. To the authors' knowledge, this property of the VMS formulation has not been dis- 
cussed/proved in the literature. We employ the standard notation used in mathematical analysis, 
for example see Reference [5]. The standard abbreviation 'a.e.' for almost everywhere is often used 
in this subsection. We now state and prove a continuous maximum-minimum principle for the VMS 
formulation. 

Theorem 3.2. Assume that the Dirichlet boundary condition is prescribed on the whole of the 
boundary (that is, = dD.), and the diffusivity tensor is assumed to be continuously differentiable. 
Let /(x) G and /(x) > almost everywhere. Let the weak solution c(x) of the variational 

multiscale formulation ()33p belong to H'^{i^) H C^{Q) H C^{0,). Then 

m_inc(x) = minc(x) 

n dn 

Proof Since c(x) G H^{n) n C^{n) n C^{n) we have 

(35) V := -nvce H^{n)nc°{Q) CV 

Define m G M and a scalar field s(x) such that 

(36) m = min c^{x.) 

(37) s(x) := max[m — c(x), 0] Vx G 

One can show that the function s(x) is piecewise C^{Q), and belongs to H^{^) n C^{d). By 
construction, we also have 

(38) s(x) > Vx G J^, and s(x) =0 Vx G r° 

Using equation (j35p (and also employing the divergence theorem) equation (j33p gets simplified 

to 

(39) (g;V.v-/) = 
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Since s(x) G H^{^) and s(x) = on F , the scalar field s(x) is a legitimate choice for (?(x). 
Substituting s(x) in the place of g(x), and noting that s(x) G H^{^1) to allow the application of 
the divergence theorem; we get 

(40) {s;vn)9n-{Vs;v)-{s-J) = 

Since = d^, s(x) = on F^, s(x) > Vx G Q, and /(x) > a.e. in Q; we conclude that 

(41) (Vs; v) = -(Vs; BVc) < 

To prove the theorem it is sufficient to show that s(x) = Vx G (which implies that c(x) > 
m Vx G il). Note that s(x) = m — c(x) unless s(x) = 0. Let 

(42) Y := {x G 17 I s(x) / 0} = {x G | s(x) > 0} 

The case Y = is trivial. We now deal with the case when Y is not empty. We first note that the 
weak derivative of s(x) is zero on 0\Y, and is — Vc on Y. This result along with equation (j4ip 
implies that 

(43) (Vs;DVs)y<0 

Since D(x) is a positive definite tensor, the above equation implies that s(x) = sq = constant 
almost everywhere in Y. Since s(x) is continuous in 0, we conclude that s(x) = sq everywhere in 
Y. Since F^ = dU C Y, and s(x) = on F^ we conclude that s(x) = on whole of Y and also on 
Y. Noting the fact that the function s(x) vanishes on the set complement of Y, we conclude that 
s(x) = on whole of Cl. Hence, we have proved the desired result. □ 

3.2. Discrete equations. The discretized finite element equations for the variational multiscale 
formulation (given by equation ([33]) ) can be written as 



(44) 







1 




\-\ 














[ h j 



where the matrices K^^ and Kpp are symmetric and positive definite, v denotes the (finite element) 
vector of nodal velocity (or auxiliary variable) degrees-of-freedom, and p denotes nodal vector of 
concentration degrees-of-freedom. Comparing the weak form (|33p and discrete equations (j44p . 
the matrices K^^, Kp^ and Kpp are obtained after the finite element discretization of the terms 
i (w; D~^v) , — (g, V • v) — i (V(/; v) and — i (Vg; DVc); respectively. Since equation ([55]) is written 
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in symmetric form, Ki,p (which comes from the term — (V • w; c) — ^ (w; Vc)) will be equal to Kp^. 
Note that, some of the terms in equation (j33|) are combined and simplified to obtain the terms 
presented in the previous line. For example, we have used the symmetry of D in obtaining the term 
— ^(w;Vc). The vectors and fp are, respectively, obtained from — (w • n;cP)pD and —{q;f) 
after the finite element discretization. 

The discrete form of equation ()34p can be written as 

(45) extremize ^-v^K^^v + p^Kp^v - \p^Kppp - v^f^ - p^ f 

V, p Z Z 

Similar to the continuous problem (that is, equations (I33p and ()34p are equivalent), the stationarity 
of the above equation (minimizing with respect to v and maximizing with respect to p) is equivalent 
to equation (I44p . By eliminating v, the Schur complement form of equation (144p can be written as 

(46) {^K p^K ^pv "l~ -^pp) P ~ -^pv^vv fv ~ fp 

Clearly, the Schur complement operator Kpi,K~^Kp^ + Kpp is symmetric and positive definite. 
The discrete variational statement of the variational multiscale mixed formulation can be posed 
solely in terms of the variable p, and takes the following form: 

(47) minimize ]-p^ [K p^K'^ K^^ + K pp) p - p"^ [K p^K'^ f ^ - f p) 

P ^ 

As mentioned earlier the variational multiscale mixed formulation does not (always) produce non- 
negative solutions for the non-negative forcing function and non-negative prescribed Dirichlet 
boundary condition. 

Remark 3.3. Unlike in the Raviart- Thomas formulation, the optimization problem (I47p is not the 
dual problem of equation (j45p . 

3.3. A non-negative formulation. A non-negative formulation based on the variational multi- 
scale formulation can be posed as the following constrained minimization problem 

(48) minimize ip^ {Kp^R-.^K^^ + Kpp) p - p^ [Kp^K'^ - f^) 

(49) subject to p ^ 

The above constrained optimization problem belongs to the class of convex quadratic programming, 
and has a unique global minimizer. 
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Remark 3.4. The variational multiscale formulation (given by equation (j33p ). in general, does 
not have the (element) local mass balance property. Specifically, one does not have the local mass 
balance property under linear equal-order interpolation for both c(x) and v(x), which is employed in 
this paper. The corresponding non-negative formulation also does not possess the local mass balance 
property. 

Remark 3.5. The non-negative method proposed in this section is also applicable for the mixed 
formulation based on the Galerkin/least- squares. As discussed in Reference [27] the variational 
multiscale and Galerkin/least- squares ( GLS) mixed formulations differ only in the definition of the 
stabilization parameter. That is, instead of the term 

i (B-^w + Vq; B (B"V + Vc)) 

which is the case for the variational multiscale formulation (see equation (|33|) ) we will have 

(D~V + Vq; t(x)D (B~V + Vc)) 

and t(x) > for the GLS mixed formulation. The discrete equations from the GLS formulation 
also takes the same form as given in equation (I44p . 

Remark 3.6. As mentioned earlier, non-negative solution is a special case of maximum-minimum 
principle. Some of the formulations presented in the literature produce non-negative solutions but 
still may violate the (general) discrete maximum-minimum principle. For example, see the non- 
negative formulation presented in Reference [30] . That is, these formulations avoid undershoots but 
may still produce overshoots. 

Though the focus of the present paper is on non-negative solutions, the proposed two non-negative 
optimization-based formulations can be easily extended to satisfy the (general) discrete maximum- 
minimum principle. To see this, let us first define the quantities Cmin niT-d Cmax to be the minimum 
and maximum values of c{x) based on the (continuous) maximum-minimum principle. Note that 
the maximum and minimum will occur on the boundary only when /(x) =0. In order to enforce 
these properties in the discrete setting modify the constraints in the corresponding optimization 
problem statements (i.e., equations l\25\i n and (j49p jgg 

(50) Cminl ^ P ^ Cmaxl 
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where 1 is a vector of ones of appropriate size. The resulting problems will still belong to qua- 
dratic programming. Hence, the proposed mathematical framework and solvers are still applicable. 
However, some interesting questions regarding the numerical performance of the solvers ( active-set 
strategy, interior point methods) need to be addressed in future work. For example, since in the case 
of general DMP we have twice the number of constraints than that in the case of the non-negative 
formulation, how large is the violation of the local mass balance in the RTO formulation because of 
the additional constraints? How much additional computational cost will be incurred because of the 
additional constraints. 

4. NUMERICAL RESULTS 

In this section we study the performance of the proposed formulations (with respect to non- 
negative solutions and local mass balance) on three canonical test problems. In our numerical 
experiments we have employed five different meshes - Delaunay, 45-degree, unstructured and well- 
centered triangular (WCT) meshes; and uniform four-node quadrilateral mesh. In two-dimensions, 
a well-centered triangulation means that all the triangles in the mesh are acute-angled (see Reference 
|31j). We will discuss more on WCT meshes in Section [4. 4[ 

The mesh layouts for the aforementioned meshes are shown in Figures [2][5l For the chosen test 
problems, the variational multiscale and Raviart-Thomas formulations in general do not satisfy 
the discrete maximum-minimum principle. We now show that, for low-order finite elements, the 
proposed two non-negative mixed formulations produce non-negative solutions for all the three test 
problems and for all the chosen computational meshes. For the variational multiscale formulation, 
we have employed equal order interpolation for the c and v fields in our numerical simulations. Note 
that (as discussed in Introduction) WCT, 45-degree, Delaunay and square meshes are sufficient to 
produce non-negative solutions for isotropic diffusion. However, these meshes may produce negative 
solutions in the case of anisotropic diffusion, which will be illustrated in this section. 

4.1. Test problem i^'S.: Anisotropic and heterogeneous medium. This test problem is taken 
from Reference [16]. The computational domain is a bi-unit square with homogeneous Dirichlet 
boundary conditions. The forcing function is taken as 



(51) / 



1 if (x,y) G [3/8,5/8]2 
otherwise 
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The diffusivity tensor is given by 

ip' + ex^ —(1 — e)xy \ 
-{l-e)xy ey^ + / 

In this paper we have taken the parameter e = 0.05. For this test problem, the numerical results for 
the concentration field c(x, y) for various meshes using the variational multiscale and corresponding 
optimization-based formulations are shown in Figure [6l The contours of the vector field v are 
shown in Figure [71 The numerical results for the concentration using the RTO and corresponding 
optimization-based formulation are presented in Figure El 

4.2. Test problem i^2: Diffusion/dispersion tensor in subsurface flows. This problem is 
taken from the groundwater modeling literature (for example, see Reference [T]). The diffusivity 
tensor is given by 

(53) B = aT||/3p + ^^Pp/3®/3 

where I denotes the second-order identity tensor, ® the standard tensor product [52], (3 the velocity 
vector, and ul and ax are respectively the longitudinal and transverse diffusivity constants. Note 
that /3 is a eigenvector of the diffusivity tensor given in equation ()53p . In this paper we have taken 
UL = 0.1 and ay = 0.01, and the velocity vector to be 

(54) /3 = e^. + ey 

where e^ and Gy are the standard unit vectors along x- and y-directions, respectively. The compu- 
tational domain is again a bi-unit square with homogeneous Dirichlet boundary conditions. The 
forcing function is same as in test problem ^1 (see equation (j51|) ). 

For the chosen velocity field (j54p (which is aligned along south-west to north-east direction) by 
rotating the current coordinate system by -|-45 degrees (i.e., in the anticlockwise direction) the 
diffusivity tensor written in the transformed coordinate system will be isotropic. Therefore, a mesh 
aligned along -|-45-degree mesh should produce non-negative solutions for the chosen velocity field 
(which is illustrated in Tables [H and . However, one will get negative solutions using a —45- 
degree mesh. For this test problem, the numerical results for the variational multiscale and RTO 
formulations and their corresponding optimization-based methods are presented in Figures [9l and 
[TOl respectively. 
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4.3. Test problem ^,^3: Non-smooth anisotropic solution. This problem is taken from Ref- 
erence [30]. The computational domain is a bi-unit square with a square hole of dimension 
[4/9,5/9] X [4/9,5/9], which is pictorially described in Figure The forcing function is taken 
as f{x) = 0. On the exterior boundary c^{x) = is prescribed, and on the interior boundary 
d'{x) = 2 is prescribed. The diffusivity tensor is given by 



In this paper we have taken ki = 1, /e2 = 100 and 9 = tt/G, which are same as the values employed 
in Reference j30j . For this test problem, the performance of the variational multiscale and RTO 
formulations and their corresponding optimization-based formulations are shown in Figure [TTl Also 
it is worth mentioning that (for this test problem) under the RTO formulation more than 45% of 
the computational domain has negative concentration, which is illustrated in Table [TJ 

4.4. Performance on a well-centered triangular mesh. In two dimensions, a well-centered 
triangulation means that each element contains its circumcenter, which is equivalent to saying that 
all elements are acute-angled triangles. A well-centered mesh in higher dimensions can be similarly 
defined [31]. Note that every WCT mesh is also a Delaunay mesh but not vice- versa. For further 
details on how to generate WCT meshes see Reference |31j . 

It is well known that well-centered triangular (WCT) meshes have some advantages in solving 
some partial differential equations as they preserve some of the underlying mathematical struc- 
ture. For example, as discussed in Introduction, a WCT mesh is sufficient to produce non-negative 
solutions for an isotropic diffusion equation. In other words, a WCT mesh respects the discrete 
maximum-minimum principle thereby preserving this key underlying mathematical property. In 
addition, a WCT mesh enables construction of a compatible discretization of a Hodge star (a ge- 
ometrical object in exterior calculus) [33]. In Reference [31] this idea has been used to construct 
a numerical method for the mixed form of the diffusion equation that is locally and globally con- 
servative, and also can exactly represent linear variation of concentration in a given computational 
domain. 

However, in this subsection we numerically show that for the RTO and variational multiscale 
formulations, even a well-centered triangular (WCT) mesh is not sufficient to produce non-negative 



(55) 
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solutions in the case of full diffusivity tensor. We consider test problem 7^1, and use the well- 
centered triangular mesh shown in Figure HI The obtained numerical results for the concentration 
using the Raviart-Thomas and variational multiscale formulations are shown in Figure [T2J As 
one can see, there are regions of negative concentration (which are indicated in white color). The 
obtained numerical results using the corresponding optimization-based formulations are also shown 
in the figure, and (as expected) we have non-negative concentration in the whole domain. 

4.5. Active set strategy and its numerical performance. The two main classes of methods 
for solving quadratic programming problems are active set strategy and interior point methods. In 
this paper we employ the active set strategy, which is very effective for small to medium sized convex 
quadratic programming. For a detailed discussion on active set strategy (including a convergence 
proof) see Luenberger and Ye |35^ Section 12.4], and for an algorithmic outline of the numerical 
method see Nocedal and Wright [24l page 462]. 

In Tables [3] and [J] we have studied the performance of active set strategy for the proposed non- 
negative formulations. We considered two cases for the initial active set. The first case is trivial, 
which is the empty set. In the second case, we have take the initial active set as the degrees- 
of- freedom that have negative concentration under the chosen formulation (either VMS or RTO). 
That is, initially one will solve a given problem using either the VMS or RTO formulation, and 
then identify the nodes (in the case of VMS) or elements (in the case of RTO) that have negative 
concentration. The initial active set is taken as those degrees-of-freedom for concentration that 
have negative values. Based on numerical experiments we found that in many problems the second 
case takes fewer iterations. However, this is not the case always, which is illustrated in Tables [3] 
and m Note that in these tables, the two choices for initial active set are denoted as 'empty set' 
and 'initial violated set.' 

4.6. Error in local mass balance under the optimization-based RTO formulation. In 
subsection 12.21 we have shown mathematically that one may have violation of local mass balance 
under the optimization-based RTO formulation. Based on the KKT optimal conditions we have 
also shown that the violation of local mass balance can occur only in the form of artificial sinks in 
some elements. These elements are those for which (^KpyV — fp)- < 0, where i denotes the element 
number. 
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In this subsection we study numerically the error in local mass balance (which is characterized 
by element sink strength). For test problems ^^1 and ^2 the total source strength is 0.0625 (which 
is equal to J^f dO). For a given element (with its boundary denoted by Fg) we calculate 
In '"^ ^ It dF, which should be negative based on the KKT conditions. (As mentioned 
earlier, in the discrete finite element setting the element source/sink strength can be obtained by 
picking the corresponding component in the KpyV — fp vector.) Contours of these element sink 
strengths are plotted using the built-in cell-centered feature in Tecplot [36] . In Figures [13] and [TJ] 
we have shown the contours of element sink strength for various computational meshes for test 
problems #1 and #2, respectively. As one can see, for these representative test problems and 
meshes the violation of local mass balance is insignificant. 

For test problem ^3 the volumetric source is zero (i.e., /(x) = in 0,). (The problem is driven 
by non-homogeneous Dirichlet boundary conditions.) Hence, we compare the element sink strength 
with the total flux along the boundary. Under the RTO formulation (that is, without optimization) 
the total (integrated) flux along the interior and exterior boundaries are —117.3852 and -1-117.3852, 
respectively. (This is not surprising as the RTO formulation has both local and global mass balance 
properties.) Under the optimization-based RTO formulation, total integrated flux along the interior 
and exterior boundaries are —117.5615 and -1-127.5694, respectively. Maximum (in magnitude) 
element sink strength is 0.7477, which is 0.636% compared to the total flux along the interior 
boundary. The total sink strength by adding the individual element volumetric (sink) strengths is 
—10.0079, which matches the difference between the fluxes along interior and exterior boundaries. 
This means that the optimization-based RTO formulation has the global mass balance property 
(but, as discussed earlier, does not possess local mass balance property). In Figure [15] we have 
shown the contours of element sink strength for test problem ^^3. 

4.7. /i- Convergence analysis. In this subsection we study the convergence of the proposed 
optimization-based methods with respect to mesh refinement. We use test problems ^^1 and ^2, 
and employ 45-degree triangular and four-node quadrilateral meshes. (As discussed earlier, we use 
-|-45-degree and —45-degree meshes for test problems ^1 and #2, respectively.) Typical four-node 
quadrilateral and 45-degree triangular meshes are shown in Figures [3] and [2]^b), respectively. We 
refine the mesh by increasing the number of nodes along each side. 
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Table 1. Performance of the RTO formulation: minimum concentration produced 
by the formulation, and percentage of elements that have negative concentrations 
(denoted as % of elements violated). 



iest problem 


1\ T 1 J- 

Mesh type 


Mm. cone. 


% of elements violated 


Problem ^1 


+45-degree 


-0.002510583 


128/648 - 


> 19.75% 




1 

Uelaunay 


-0.000011158 


67/1800 - 


3.72% 




Well-centered 


-0.000824908 


45/336 - 


> 13.39% 


Problem ^2 


-|-45-degree 


0.000000000 


0/648 - 


0.00% 




-45-degree 


-0.006674991 


216/648 - 


> 33.33% 




Delaunay 


-0.000468342 


71/1800 - 


3.94% 




Well-centered 


-0.000616864 


42/336 - 


> 12.50% 


Problem ^^3 


Mesh in Figure [5] 


-0.081453689 


848/1868 - 


> 45.40% 



Table 2. Performance of the variational multiscale formulation: minimum concen- 
tration produced by the formulation, and percentage of nodes that have negative 
concentrations (denoted as % of nodes violated). 



Test problem 


Mesh type 


Min. cone. 


% of nodes violated 


Problem #1 


-|-45-degree 


-0.000986744 


30/361 - 


8.31% 




Delaunay 


-0.000000531 


11/961 - 


1.14% 




Well-centered 


-0.000056615 


4/191 - 


2.09% 




Four-node quadrilateral 


-0.000000901 


7/361 - 


1.93% 


Problem ^2 


-|-45-degree 


0.000000000 


0/361 - 


0.00% 




-45-degree 


-0.000402642 


40/361 - 


> 11.08% 




Delaunay 


-0.000000941 


12/961 - 


1.25% 




Well-centered 


-0.000007018 


5/191 - 


2.62% 




Four-node quadrilateral 


-0.000000155 


6/361 - 


1.67% 


Problem #3 


Mesh in Figure [5] 


-0.004613415 


264/998 - 


> 26.45% 



In Tables [5] and [6] we show the variation of minimum concentration and percentage of nodes 
that produce negative solutions with respect to mesh refinement under the VMS formulation for 
45-degree triangular and four-node quadrilateral meshes. In Table [7] we have shown the variation 
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Table 3. Performance of the active-set strategy for the RTO formulation. We have 
employed the D2-RT0 problem given by equation ([25]) . 



iest problem 


Mesh type 


Active set iterations 
(empty set, mitial violated setj 


Problem #1 


-|-45-degree 


(lU5,43j 




Delaunay 


(25,53) 




Well-centered 


(31,20) 


Problem #2 


-|-45-degree 


N/A 




-45-degree 


(185,57) 




Delaunay 


(29,51) 




Well-centered 


(32,15) 


Problem ^^3 


Mesh in Figure [5] 


(658,436) 



Performance of the active-set strategy for the variational multiscale for- 
We have employed the non-negative formulation given by equations (j48p 



Table 4. 
mulation. 
and dMl). 



Test problem 


Mesh type 


Active set iterations 
(empty set, initial violated set) 


Problem ^\ 


+45-degree 


(18,15) 




Delaunay 


(8,5) 




Well-centered 


(4,2) 




Four-node quadrilateral 


(7,3) 


Problem #2 


-|-45-degree 


N/A 




-45-degree 


(25,17) 




Delaunay 


(6,9) 




Well-centered 


(5,2) 




Four-node quadrilateral 


(7,1) 


Problem 


Mesh in Figure [5] 


(72,240) 
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of minimum concentration and percentage of elements that produced non-negative solutions with 
respect to mesh refinement under the RTO formulation. (Note that in the RTO formulation, the 
concentration, c(x), is a constant over each element.) As expected, the percentage of non-negative 
nodes and minimum concentration decrease with respect to mesh refinement, but still the persis- 
tent non-negative values and spatial extent of the violation prohibits the use of VMS and RTO 
formulations in simulations for which transport is coupled with reactions. 

In Figures [TBI and [T71 we have compared the number of iterations taken by the active set strategy 
with respect to mesh refinement for the proposed optimization-based methods. As one can see 
from these figures, the number of iterations stabilized with respect to mesh refinement (that is, the 
strategy takes almost the same number of iterations as the mesh is refined). The optimization- 
based VMS method takes fewer active-set strategy iterations compared to the iterations taken by 
the optimization-based RTO formulation. 

Figures [18] and [191 respectively, compare the CPU time taken by the optimization-based VMS 
and RTO methods with respect to mesh refinement. On the y-axis we plot the ratio between the 
additional CPU time taken by the active-set strategy (or the optimization solver to obtain non- 
negative solution) and the CPU time taken by the corresponding underlying mixed formulation 
(either VMS or RTO formulation). For the optimization-based VMS method, the additional cost 
to obtain non-negative solution using the four-node quadrilateral mesh is only a fraction of the 
computational cost of the VMS formulation. For the three-node triangular mesh, the additional 
cost to obtain the non-negative solution is nearly twice the cost of the VMS formulation. The 
optimization-based RTO method takes relatively more additional CPU time to obtain non-negative 
solution, and the ratio between the additional CPU time and the CPU time taken by the RTO 
formulation is nearly 5 for test problem ^1 and 20 for test problem ^2. This should not be 
surprising as in the RTO formulation the primary variable (in our case, the concentration) is poorly 
approximated by its piecewise constant representation over each element. (Note that, in the VMS 
formulation the primary variable is continuous, that is, piecewise linear and continuous across 
elements.) 

In Figure [20] we compare the error in the local mass balance with respect to mesh refinement 
for the RTO formulation. Since, in the non-negative version of the RTO formulation we always 
have artificial sinks (that is, the error in the local mass balance will always be negative), we have 
taken the negative of the error before taking the logarithm. We have plotted both the maximum 
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(artificial) element sink strength (or error in local mass balance in an element) and also the total 
sink strength by summing the contribution from all elements. From Figure [20] one can see that, for 
the chosen problems, the error in the local mass balance decreases exponentially with respect to the 
element size. 
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Table 5. Performance of the variational multiscale formulation using three- node 
triangular element with respect to mesh refinement. We have employed +45-degree 
and — 45-degree meshes for test problems #1 and #2, respectively. 



Test problem 


# of nodes per side 


Mm. L-onc. 


% nodes violated 


Problem ^1 


1 n 
iU 


-O.iyilj-UUo 


9/100 


9% 




19 


-9.87E-004 


30/3d1 


010/ 

8.3170 




2o 


1 ni T7^ r\r\A 


54/784 ^ 


6.89% 




o / 


Q AIT? nnfi 

-o.4iIl/-UUD 


64/1369 ^ 


4.67% 




40 


-LUOHj-UUD 


71/2116 ^ 


3.36% 




00 


-Z.L)511j-UU ( 


71/3025 ^ 


2.35% 




04 


V QQT7 nno 


75/4096 ^ 


1.83% 






-o.4di1j-UUo 


77/5329 ^ 


1.44% 


Problem ^2 




1 OQT7 nno 


6/100 


^ 6% 




19 


-4.(Jo11j-UU4 


40/361 ^ 11.08% 




28 


-4.72E-005 


66/784 ^ 


8.42% 




37 


-5.41E-006 


66/1369 ^ 


4.82% 




46 


-9.57E-007 


66/2116 ^ 


3.12% 




55 


-2.40E-007 


66/3025 ^ 


2.18% 




64 


-8.50E-008 


66/4096 ^ 


1.61% 




73 


-3.45E-008 


66/5329 ^ 


1.24% 



5. CONCLUSIONS 

Tensorial diffusion problems arise in a variety of important engineering and scientific applica- 
tions. Although the continuous problem satisfies a maximum-minimum principle, most numerical 
approximations fail to satisfy this principle in a discrete sense on arbitrary meshes. For some appli- 
cations, violation of the discrete maximum-minimum principle can be problematic due to physically 
meaningless negative values of the dependent variable. 

In this paper, we proposed two non-negative low-order mixed finite element formulations for 
the tensorial diffusion equation. (That is, the proposed formulations provide non-negative numer- 
ical solutions for linear, bilinear and trilinear finite elements.) This is achieved by rewriting the 

26 



Table 6. Performance of the variational multiscale formulation using four-node 
quadrilateral element with respect to mesh refinement. 



Test problem 


jj^ of nodes per side 


Min. cone. 


% of nodes violated 


Problem ^1 


ii 


-i.i0rj-UU4 


7/121 - 


5.7970 






-z.UzJij-UU / 


8/441 - 


1.8170 




Q1 




o 1 (\r' 1 

o/ybi - 


n o oO/ 

U.8c)7o 




A 1 

41 


-i.Uiil;-UUy 


9/1681 - 


0.54% 




51 


-2.ozr!j-UiU 


9/2601 - 


0.35% 


Problem #2 


11 


-1.15E-004 


6/121 - 


4.96% 




21 


-2.02E-007 


6/441 - 


1.36% 




31 


-7.17E-009 


6/961 - 


0.62% 




41 


-l.OlE-009 


6/1681 - 


0.36% 




51 


-2.32E-010 


6/2601 - 


0.23% 



Table 7. Performance of the RTO formulation with respect to mesh refinement. 
We have used -|-45-degree and — 45-degree meshes for test problems #1 and ^^2, 
respectively. 



Test problem 


^ of nodes per side 


Min. cone. 


% of elements 


violated 


Problem ^1 


10 


-0.029349220 


48/162 - 


> 29.63% 




19 


-0.002510583 


128/648 


> 19.75% 




28 


-0.000065471 


158/1458 


> 10.84% 




37 


-0.000027412 


194/2592 - 


7.48% 




46 


-0.000012794 


224/4050 - 


5.53% 




55 


-0.000006529 


248/5832 - 


4.25% 


Problem #2 


10 


-0.039682770 


64/162 - 


> 39.51% 




19 


-0.006674991 


216/648 


> 33.33% 




28 


-0.001901262 


384/1458 


> 26.34% 




37 


-0.000756689 


504/2592 - 


> 19.44% 




46 


-0.000337093 


552/4050 


> 13.63% 




55 


-0.000140215 


576/5832 - 


9.88% 
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formulations as constrained optimization problems. In both the cases, the problem belongs to con- 
vex quadratic programming, which can be effectively solved using existing numerical optimization 
solvers (e.g., active set strategy and interior point methods). 

One of the formulations is based on the variational multiscale formulation, and the other is based 
on the lowest-order Raviart-Thomas spaces (that is, the RTO element). We have demonstrated 
that the variational multiscale formulation satisfies a continuous maximum-minimum principle. In 
the case of non-negative formulation based on Raviart-Thomas spaces, two different optimization 
problems are presented - the primal and dual problems. From the optimization theory it has 
been inferred that these two problems are equivalent. In addition, from the Karush-Kuhn- Tucker 
optimality conditions it is inferred that one may have violation of (element) local mass balance in 
those elements that have zero concentration. These violations of local mass balance are in some 
limited part of the domain, and the violations are small for the test cases studies here. 

We have studied the convergence properties of the proposed optimization-based methods. Through 
numerical experiments we have shown that the error in local mass balance under the optimization- 
based RTO method decreases exponentially with respect to mesh refinement. We have also studied 
the performance of active-set strategy method for solving the resulting convex quadratic program- 
ming problems. The performance of the proposed non-negative formulations is illustrated on three 
representative problems, and the formulations performed well. 

One note worthy feature is that existing solvers based on the variational multiscale and RTO 
formulations can be easily extended to implement the proposed optimization-based non-negative 
formulations. Designing a non-negative (stabilized) mixed formulation that also possesses local 
mass balance property is part of our future work. 

APPENDIX: SOME CLASSICAL RESULTS ON DISCRETE MAXIMUM-MINIMUM 

PRINCIPLE 

One of the early works on DMP dating from the 1960s is by Varga |371I38|. and was presented in 
the context of finite difference schemes. To the authors' knowledge, the initial work on DMP in the 
context of the finite element method was done by Ciarlet and Raviart [Sj. (Another relevant work 
by Ciarlet is [39], which was in the context of finite difference operators.) Reference [8] addressed 
linear simplicial finite elements, and considered the classical Galerkin single-field formulation for 
the Poisson and Helmholtz equations. 
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As mentioned earlier, classical finite element formulations do not satisfy the DMP on general 
meshes for full diffusivity tensor. The formulations which satisfy DMP impose severe restrictions 
on both meshes and coefficients of the diffusivity tensor. We now outline some of the classical 
results that are available on DMP in the context of the finite element method for a scalar diffusion 
equation. 

• For linear simplicial elements, Ciarlet and Raviart [8] have shown that non-obtuseness is 
sufficient to satisfy DMP. 

• Christie and Hall [ID] have presented sufficient conditions for bilinear finite elements to 
satisfy the DMP under a homogeneous forcing function. The results can be summarized 
as follows. For a non-uniform rectangular mesh (see Figure [T]), the DMP will be satisfied 
provided 

(56) /11/12 ^ 2™^"^ ('^i' ^2) kik2 < — max(/i^, /ig) 

This implies that a uniform rectangular mesh (that is, hi = /i2 and /ci = /C2) will satisfy the 
DMP provided 

(57) -^ki <hi<V2ki 

This further implies that a mesh with squares (that is, hi = h2 = ki = /C2) will always 
satisfy the DMP for the case of a scalar diffusion equation. 

• Vanselow [H] has shown that a Delaunay triangulation along with an additional condition 
on boundary nodes are sufficient for the DMP under the classical single-field Galerkin 
formulation. We now outline how this additional condition looks for a convex domain. 
To this end, let Pi and P2 be two neighboring nodes on the boundary. Let us denote 
their spatial coordinates as xi and X2, and define x := (xi +X2) /2. Then the additional 
condition can be written as 

(58) ||xj — x||2 < ||xq — x||2 for all nodes Q in the triangulation with Q ^ Pi, i = 1,2 

where xq denotes the spatial coordinates of the node Q. A similar condition is required for 
a non-convex domain. In addition, it has also been shown that under some weak additional 
assumptions on the triangulation, these conditions are necessary [HI Section 4]. 
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Figure 1 . Non-uniform rectangular mesh where hi and /i2 denote the length of hor- 
izontal edges of two neighboring elements. Similarly, ki and k2 are for corresponding 
vertical sides. 
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(a) Delaunay mesh 




0.2 0.4 0.6 0.8 1 

X 

(b) +45-degree mesh 



Figure 2. Typical triangular meshes [Delaunay (top) and +45-degree (bottom)] 
that are used for test problems 1-3 are shown in the figure. In addition, — 45-degree 
mesh is also used in numerical simulation^^^ which is similar to the -|-45-degree except 
that the diagonals run along south-east to north-west direction. 
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Figure 5. Pictorial description of test problem ^2>: Computational domain is the 
bi-unit square with a square hole [4/9,5/9] x [4/9,5/9]. On the exterior boundary 
cP{x) = is prescribed. On the interior boundary (?{x) = 2 is prescribed. The 
computational domain is triangulated using Gmsh |42j . 
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Figure 6. Test problem ^1: Performance of the variational multiscale (left) and 
corresponding optimization-based (right) formulations. In the numerical simulations 
45-degree (top), Delaunay (middle) and four- node quadrilateral (bottom) meshes are 
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employed. The regions that have negative concentration are indicated in white color. 
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Figure 7. Test problem Contours of the components of the vector field v ob- 
tained using the variational multiscale (left) and corresponding optimization-based 
(right) formulations. The top figures show the x-component (that is, f^), and the 
bottom figures the y-component of v. Four-node quadrilateral mesh is used in the 
numerical simulation. 
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Figure 8. Test problem #1: Performance of the RTO triangular (left) and cor- 
responding optimization-based (right) formulations. In the numerical simulations 
-|-45-degree (top) and Delaunay (bottom) meshes are employed. The regions that 
have negative concentration are indicated in white color. 
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Figure 9. Test problem Performance of the variational multiscale (left) and 
corresponding optimization-based (right) formulations. In the numerical simulations 
Delaunay (top), — 45-degree (middle) and four-node quadrilateral (bottom) meshes 

40 

are employed. The regions that have negative concentration are marked as white. 
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Figure 10. Test problem ^2: Performance of the RTO triangular (left) and cor- 
responding optimization-based (right) formulations. In the numerical simulations 
Delaunay (top) and — 45-degree (bottom) meshes are employed. The regions that 
have negative concentration are indicated in white color. 
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Figure 11. Test problem jj^?,: Performance of the variational multiscale (top) and 
RTO (bottom) formulations. The figures on the right show the performance of cor- 
responding optimization-based non-negative formulations. As expected, the varia- 
tional multiscale and RTO formulations produced negative solutions in significant 
portions of the domain, which are denoted by the white color. On the other hand, 
the proposed optimization-based formulations produced desirable non-negative solu- 
tions. The computational mesh that is used in these numerical simulations is shown 
in Figure [3 
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Figure 12. Performance on a well-centered triangular (WCT) mesh: RTO (top) 
and variational multiscale (bottom) formulations on test problem ^1 and using 
the WCT mesh shown in Figure [H The left figures illustrates that, in the case 
of tensorial diffusivity coefficient, even a WCT mesh is not sufficient for the RTO 
and variational multiscale formulations to satisfy the discrete maximum-minimum 
principle. The right figures shows that (as expected) the proposed non-negative 
optimization-based formulation outlined in Sections [2] and [3] produces non-negative 
solutions. 
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Figure 13. Test problem ^1: Contours of mass balance error under the 
optimization-based RTO formulation. In the numerical simulations +45-degree 
(top), Delaunay (middle) and WCT (bottom) meshes are employed. 
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Figure 14. Test problem ^2: Contours of mass balance error under the 
optimization-based RTO formulation. In the numerical simulations -45-degree (top), 
Delaunay (middle) and WCT (bottom) meshes are employed. 
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Figure 15. Test problem ^^3: Contours of mass balance error under the 
optimization-based RTO formulation. 
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Figure 16. This figure compares number of iterations taken by the active set 
strategy for the VMS formulation with respect to mesh refinement. We have em- 
ployed four-node quadrilateral (top figure) and three-node triangular (bottom figure) 
meshes. Note that, in the case of three- node triangular meshes, we have used +45- 
degree and — 45-degree meshes for test problems #1 and #2, respectively. For both 
the test problems two different sets are used as an initial guess in the active set 
strategy. 
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Figure 17. This figure compares number of iterations taken by the active set strat- 
egy for the RTO formulation with respect to mesh refinement. We have employed 
+45-degree and — 45-degree triangular meshes for test problems #1 and #2, respec- 
tively. For both the test problems two different sets are used as an initial guess in 
the active set strategy. 
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Figure 18. This figure compares the computational effort of the optimization-based 
VMS method with respect to mesh refinement. We employed four- node quadrilateral 
(top figure) and 45-degree triangular (bottom figure) meshes. On the y-axis we have 
the ratio of the additional CPU time taken by the optimization-based VMS method 
to the CPU time taken by the VMS method (which produces negative solutions). 
We considered test problem ^1 and ^2. Note that each node has three degrees-of- 
freedom. 
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Figure 19. This figure compares the computational effort of the optimization-based 
RTO method with respect to mesh refinement. On the y-axis we have the ratio of the 
additional CPU time taken by the optimization-based RTO method to the CPU time 
taken by the RTO method (which produced negative solutions). We employed +45- 
degree and — 45-degree triangular meshes for test problems #1 and #2, respectively. 
Note that each node has three degrees-of-freedom. 
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Figure 20. This figure compares the error in the local mass balance with respect 
to mesh refinement using the RTO formulation. We have used +45-degree and 
— 45-degree triangular meshes for test problems #1 and #2, respectively. We have 
plotted both the maximum element sink strength, and the total error by summing 
the artificial sink strengths in all elements. From the above figure one can see that 
both these errors in the local mass balance decrease exponentially with respect to 
the element size. 



51 



